Transport in a highly asymmetric binary fluid mixture 



Sorin BasteaB 

Lawrence Livermore National Laboratory, P.O. BOX 808, Livermore, CA 94550 

. - We present molecular dynamics calculations of the thermal conductivity and viscosities of a model 

colloidal suspension with colloidal particles roughly one order of magnitude larger than the suspend- 
^— »' ' ing liquid molecules. The results are compared with estimates based on the Enskog transport theory 

^ ' and effective medium theories (EMT) for thermal and viscous transport. We find, in particular, 

that EMT remains well applicable for predicting both the shear viscosity and thermal conductivity 
O ■ of such suspensions when the colloidal particles have a "typical" mass, i.e. much larger than the 

D ' liquid molecules. Very light colloidal particles on the other hand yield higher thermal conductivi- 

, ties, in disagreement with EMT. We also discuss the consequences of these results to some proposed 

mechanisms for thermal conduction in nanocoUoidal suspensions. 

PACS numbers: 66.60.+a, 66.20.+d, 82.70.Dd, 83.80.Hj 

D ' Liquid suspensions of solid particles (colloids) are widely encountered in biology, industry and many natural pro- 
cesses. In addition to their relevance for numerous practical applications they have emerged as a useful paradigm 
for the study of phase transitions, from crystal nucleation and growth to gelation Colloids have not only inter- 
esting thermodynamic properties, but remarkable rheological properties as well i and very complex flow behavior 
When the suspended particles are only one to two orders of magnitude larger than typical liquid molecules, i.e. 
in the nanometer domain, colloids may exhibit entirely new properties 0, Q that are expected to have important 
I technological consequences. A theoretically important and practically relevant class of colloids consists of suspensions 
of spherical colloidal particles with interactions dominated by excluded volume effects. In the following we employ 
molecular dynamics simulations to calculate the thermal conductivity and viscosities (shear and bulk) of fairly dilute 
■ O ■ colloidal suspensions modeled as mixtures of strongly asymmetric spherical particles interacting through short range 
repulsive potentials. We discuss the results in the light of theoretical estimates based on microscopic and macroscopic 
pictures of the system. Despite the simplicity of the model the conclusions should provide, inter alia, some guidance on 
the expected transport properties of dilute suspensions of nanosized particles, particularly the thermal conductivity, 
which has been the subject of some speculation 0, H, 0|- 

The model that we study consists of two types of particles, 1 - solvent and 2 - colloid, with masses mi = m and 
TO2 = mc. The interaction potentials between the particles are based on the inverse-12, "soft sphere" potential, whose 
properties have been well studied [1, Q 

g : = ^ (7) (1) 

, and which we truncate and shift at r/cr = 2; we also define u{r) = 00 for r < 0. The interactions are: 

; uii(r)=u(r) (2) 

S; ui2ir) ^ uir - R,) (3) 

U22{r) = u{r - 2Rc) (4) 

ch : 

Q , Similar interactions, that take into account the "size" of the colloidal particles by introducing a "hard core radius" 
O ^ Rc, have been employed before to model suspensions [l^, |Tl[. For temperatures ksT ~ e the effective diameters 
^ I corresponding to the above interactions should be well approximated by ai = a, <ti2 = Rc + <7, and (72 = = 2Rc + a, 

■ and satisfy (T12 — (cti -I- tT2)/2. In the following we will therefore quote as relevant quantities the "diameter" ratio of 

■ the colloid and solvent particles, CTc/ct, and the "volume fractions" of the colloidal particles, (j)c = Trnc(Tc^/6, and the 
^ . solvent, = Tma'^/6; Uc and n are the corresponding number densities, Uc = Nc/V, n = N/V. All the simulations 

■ - - ' presented here were performed in the microcanonical (NVE) ensemble with the average temperature set to ksT — e. 

The systems (mixtures) that we studied are summarized in Table 1, and correspond to two "diameter" ratios, 
each with two mass ratios. Since for a realistic colloidal particle the ratio of its mass to that of a fluid molecule 
is nic/m (ac/cr)^, our perhaps most practically relevant results correspond to rric/m — 1000. However, we also 
analyzed the effect of a much smaller mass ratio, vrLc/m — 1. The volume fractions ipc of the colloidal particles have 
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been chosen low enough so that the system is rather dilute, but high enough so that a reasonable number of colloidal 
particles can be simulated without the need for a prohibitively large number of solvent particles. (Nevertheless, N is 
rather large, of the order 10^.) The pairs cj) and 0c for the two different diameter ratios have been chosen to yield the 
same system pressure poj corresponding to a pure solvent at UQa^ — 0.8 {(fio — 0.419). Incidentally, we have found 
that this can be accomplished with very good precision (better than 1%) by using the scaling relation: 

p(g, 0c) ^ PMCSl{<7, gc, 0, 0e) 

Po(ct, 00) Pcs{(^Ao) 

where p is the system pressu re, p bmcst, the Boublfk-Mansoori-Carnahan-Starling-Leland equation of state pressure 
of a hard-sphere mixture (l^ . Il3l | and pes the Carnahan-Starling equation of state pressure of the hard-sphere liquid 
. The choice of a common pressure allows an unambiguous comparison of the thermal conductivity and viscosities 
of the suspension with that of the reference system (pure solvent at pressure po), is implicitly assumed by theories 
relying on macroscopic scale arguments (see below), and should also correspond to the usual experimental situations. 

The calculation of the viscosities (shear and bulk) and thermal conductivity can be done in molecular dynamics 
simulations using the Green-Kubo relations, which express these linear transport coefficients as time integrals of auto- 
correlation functions of microscopic currents [isl . [l^ . . This formalism yields unambiguous definitions for the shear 
and bulk viscosities, applicable to both single fluids and mixtures: 

1 /■* ^ 

- ^77^ / ( E ^o.0yo.fi{T))dT (6a) 
m = HTT^ E Ka(0) -P\)dr (6b) 

^ "'0 a,/3=l 

where a is the microscopic stress tensor: 

o'q/3(t) = E [^iVia{'r)vii3{'r) -|- (T)rj^ (t)] (7) 

i 

(a, (3 = x,y, z), p is the pressure, and the viscosities rj and ^ are given by the i — > oo limits of the above relations. 

The treatment of thermal transport in mixtures on the other hand is more complicated due to the coupling of 
energy and mass transport Since this is an important but many times confusing issue we discuss it briefly below 
for the present case of a binary mixture (see also the discussion in [Tsj). The hydrodynamic equations for a binary 
mixture express species conservation, as well as momentum and entropy transport: 



V-(paVa)=0 (8a) 



dpa 

dt 
dv 

p— -I- pv Vv = -V • P (8b) 

p-^ + pv • Vs = -V • Js + e (8c) 

In the above equations pa and Va (a = 1, 2) are the (position and time dependent) mass densities and flow velocities of 
the two species, respectively; p is the total mass density, p = pi +p2; v is the center of mass ( "barycentric" ) velocity, 
V = p~^(piVi -I- P2V2); P is the stress tensor, Pai3 = pSafi — P'afi, with p hydrostatic pressure and P/ viscous stress 
tensor, P'a0 = [ri{dva/ dxp + dvp/dxa) + (C ~ 2?7/3)V • v5q,^]; s is the entropy density, the entropy current and 
the entropy production. The entropy current is expressed in terms of heat - J^, and mass diffusion - Jq, currents: 



3s = — [Jq - (miJi + P232)] (9a) 



1 

J, = - (pev + P • v) (9b) 

Ja=Pa(Va-v) (9c) 



, with Pa chemical potential (per unit mass), e total specific energy and Jg the corresponding energy current. In 
the framework of non-equilibrium thermodynamics [13] the heat and mass currents (denoted as the set {J5}) are 
connected to thermodynamic forces {X5} by heat-mass linear transport coefficients {La^}, ((5,7 — 1,2,(7): 

ls=Y,Ls^^^ (10) 
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The entropy production O contains independent contributions from "vectorial" phenomena (heat and mass transport) 
- Qy and "tensorial" ones (momentum transport) - Ot, O ~ G„ + Qf. Both contributions assume the Onsager form, 
i.e. for heat-mass processes: 

s 

while Qt — (l/r)P/ : Vv. The Onsager reciprocity relations, Ls^ — L^g, along with Ji + J2 = from the definition 
of the diffusion currents leave 3 independent heat-mass transport coefficients for a binary mixture, {L12, Liq, Lqq}, 
with the currents now written as 

Ji = -Li2(Xi -X2)+Li,X, (12a) 

J, =Llg(Xl -X2)+LggX, (12b) 

Three distinct sets of currents and thermodynamic forces have been discussed in detail [TtI . [isj , each with different 
transport coefficients. The currents (and forces) of different sets are connected by linear transformations under which 
is invariant and preserves its Onsager form. This leads to well defined and useful relations between the coefficients 
p^ . We would like to point out that one such relation can be deduced without considering in detail the particular 
definitions of currents and forces. We simply note that all sets use the same, physically intuitive, heat driving force, 
Xg = — VT/T, as well as diffusion currents Ja given by Eq. 9c. Since the phenomenological definition of the thermal 
conductivity A is based on the observation that in the absence of diffusion, i.e. Ji = 0, the heat current should reduce 
to its canonical form, = — AVT [l^, this yields 

+ (13) 



It is worth noting that, as opposed to Liq and Lqq, X does not depend on the chosen set of currents and forces and 
moreover, it assumes the above form for all the sets. 

For our calculations we adopt the "mainstream" choice for forces (and currents), 

X„ = -TV {^la/T) (14a) 

X, - -^VT (14b) 

(a — 1,2) 18], but this selection is not in fact arbitrary. As first discussed by Erpenbeck [l^, the "mainstream" 
set is preferable for molecular dynamics calculations since its corresponding microscopic currents only depend on 
microscopic quantities easily available in simulations. The other choices on the other hand require the knowledge of 
thermodynamic quantities such as chemical potentials or partial enthalpies which are difficult to calculate with any 
accuracy (see also below). 

The microscopic currents for the "mainstream" set are 

= E'^*['^^('^) ~ Vcm{t)] (15) 
i(a) 

+ ^ E E [r^(^) - ^^■(^)] ^*(^) • (^) - H^CMir) (16) 

where H and m are the enthalpy and center of mass velocity of the system, respectively. Since vcm is set to zero 
in the simulations, H does not enter in fact the calculations. The Green-Kubo relations for the heat-mass coefficients 
are: 



L,S) = ^^7^ (j«(0) • j«(r))dr (17a) 



Li,(t) = T^TTTTTF / (ji(0) ■hir))dT (17b) 



3VkBT Jo 
1 

SVkBT 



Li2{t) - TTTrt^ / (ji(0) ■h{r))dT (17c) 
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and {Li2, Lig, Lqq} correspond to the t — > oo limits of the above relations. 

It has been sometimes remarked [20J that the "mainstream" set does not allow a proper calculation of the thermal 
conductivity since A as defined by Eq. [13] may result from the subtraction of two large quantities leading to significant 
errors. To avoid this perceived problem these authors used a different set of currents, with microscopic heat current 

i'q = iq - C'^lh + h2h) (18) 

where hi^2 are the partial specific enthalpies. This was expected to shift most of the thermal conductivity contributions 
to the first term of Eq. 1131 and therefore result in a more appropriate definition. We would like to remark that, if we 
define a time-dependent thermal conductivity 



Lqq{t) 



Llq{t)Lqlit) 



(19) 



which satisfies A = limt^ao^{t) , it is easy to show that A(i), similarly with A, is also invariant under such a change of 
currents. Consequently, to the extent that the thermal conductivity is defined as usual from the long-time "plateau" 
of A(i), using the heat current given by Eq. [TH]does not offer any real advantage. 

The molecular dynamics calculation of the transport coefficients relies on Eqs. [n]and[T71 whose integrands are easily 
calculated during simulations. We have performed such calculations for the systems described in Table 1 and also the 
reference (pure solvent) system. The units for viscosity and thermal conductivity have been chosen (mA:BT)^/^/(T^ 
and (fc^T/m)^/^/(T^, respectively; the time unit is to = a{m/kBTY^'^ ■ In these units we find that the reference system 
has viscosities ryo = 1.11, Co = 0.21 and thermal conductivity Aq = 4.87. To provide an intuitive connection to the 
often studied hard sphere system we also estimate for the reference system a mean free time between "collisions" of 
r ~ 0.035. 

Some of the autocorrelation functions (integrands) for the viscosities and thermal conductivity (Eqs. 6 and 17) are 
shown in Figs. 1 and 2 (normalized by their t = values) for the reference system and the colloidal system with 
the largest size ratio, Uc/'J — 15. One interesting feature of the shear and bulk viscosity integrands is that they are 
largely independent of the mass ratio mc/m-, even when it varies by three orders of magnitude. This feature also 
extends to the transport coefRcients themselves; as shown in Fig. 6 the shear viscosity of the colloidal system with 
mc/m = 1000 appears to be largely similar with that of the mc/m = 1 system. Despite the fairly small colloidal 
"volume fraction" the time-dependent shear viscosity ij{t) exhibits for both (Tc/cr = 15 systems a pronounced early 
times peak, corresponding to a significant viscoelastic response [3]. This effect appears much reduced when the 
diameter ratio is only slightly smaller, adu = 10 - Fig. 6. The same behavior is also observed for the bulk viscosity 
- Fig. 7. 

The heat-mass autocorrelations - Figs. 3-5, behave qualitatively very different from the viscosity ones as a function 
of the mass of the colloidal particles. We note for example that they exhibit strong oscillations for light colloidal 
particles, i.e. for rric = m, and a much smoother character for heavy ones, i.e. for mc/m — 1000. The thermal 
conductivity \{t) itself reflects these differences both at early times and as i — > oo. 

In the following we would like to compare the MD results for 77, C and A with available theories for transport in 
suspensions. The prediction of the transport properties of the present model colloidal suspension can proceed in 
principle along two different paths. The flrst path views the suspending liquid as a structureless matrix (continuum) 
and the colloidal particles as "impurities" (or "dispersed phase") with well defined properties distinct from those of 
the matrix. Then, by evaluating the response of the system to small, macroscopically applied fields, e.g. large scale 
temperature gradients or imposed shear fiows, the equivalent, effective transport properties of the system can be 
determined (see, for example, Refs. [21I, H^l). This method has a long history and is commonly known as effective 
medium theory (EMT). It has been successfully applied to both liquids and solids containing "impurities" which 
are large compared to any inherent matrix structure and sufficiently far away from each other , i.e. dilute. The 
predictions of these theories typically depend only on the volume fraction (j) occupied by the dispersed phase [23j . 
as well as "matrix" and "dispersed phase" properties. For dilute enough systems the (f) dependence is with a good 
approximation linear. When applied to the present case, where the colloidal particles considered are both "solid" and 
thermally " insulating" , the transport coefficients will therefore be functions of the liquid " matrix" properties alone, 
arguably at the same pressure and temperature. Such relations have been in use for more than a century, and yield 
for the suspension viscosity [l^, i21„ 2^ , 

Veff=Ml + l^) (20) 

(Einstein's result for the viscosity of a dilute suspension), while for the thermal conductivity [2^ 

Ae//=Ao(l-^0) (21) 
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For the bulk viscosity the only similar result that appears to be available [25j suggests that there is no contribution to 
the effective bulk viscosity to first order in cf). Considering the above relations, the application of EMT to the present 
system would therefore seem to be straightforward and require only the value of the volume fraction (j) occupied by 
the impurities, i.e. colloidal particles. Unfortunately this value is rather ambiguous when the modeling is done at the 
microscopic level, particularly when the size of the impurities is comparable with the size of the fluid particles, as is 
the case here. The problem is that EMT interprets 1 — </) as the volume fraction occupied by the fluid matrix, which 
is itself equivocal. This volume fraction may be defined as either 1 — (j)c, or perhaps better for small enough (/)c, as 
1 — where (pfc = 0c(l + cr/ac)^ ^IQj : 0/c accounts in the usual way for the actual exclusionary volume around a 
colloidal particle, where the presence of a solvent molecule is forbidden (in hard spheres language). Any reasonable 
definition of the fluid accesible volume fraction would arguably be bounded by these two values, 1 — </> and 1 — (j)/,,, 
and we therefore quote effective medium theory predictions corresponding to both of them. Although the difference 
between 0c and 0/c is not completely negligible, this does not affect our conclusions. 

The second, conceptually different treatment of transport in suspensions regards the system as a binary fluid 
mixture, i.e. it considers its microscopic, particle character and its detailed interparticle interactions. While no fully 
microscopic theory for the transport coefficients of either simple fluids or mixtures is available, the Enskog theory 
(ET) for the hard sphere fluid has proved to be successful in a significant thermodynamic domain [26], and with 
suitable modifications has been shown to be applicable to other relevant simple fluids [13, . The corresponding 
theory for hard sphere mixtures has also been rigorously derived poj , and tested using MD simulations for a number 
of relevant cases [18, 30, 31]. Although the binary mixture studied here is modeled by "soft sphere"-based potentials 
and not hard spheres, we employ a simple scaling relation to estimate the relative values of the mixture transport 
coefficients with respect to those of the reference system at the same pressure and temperature, i.e. 

H(a, (Tc, 0, (f>c, m, mc)/Ho(cr, (f>Q,m) = E,et{o; ctc, 0, (f>c, m, mc)/H_ETo(o-, (t>o,m), (22) 

where H stands for 77, ^ or A, and denotes the reference (pure solvent) system. The Enskog theory hard sphere 
results (denoted above by ETq) are well known [2^, while the Enskog mixture theory relations can be found in [29| : 
they are too complicated to be meaningfully quoted here. We used the second Enskog approximation [2^ and tested 
our numerical implementation against the values quoted in [sT]. 

We now proceed to compare the MD simulation results for shear and bulk viscosities and thermal conductivity with 
the predictions of these two theories - see Table 2. (We note first that the quoted MD values carry statistical error bars 
of approximately 10-15%, which should also encompass small deviations due to the neglect of long time tails [sol-lsH.! 
The MD calculated shear viscosity is fairly well reproduced by the effective medium theory Eq. [201 for both large and 
small colloidal masses. The Enskog predictions on the other hand appear to deviate significantly from the MD values, 
particularly for the larger colloidal particles. Surprisingly, the opposite seems to hold for the bulk viscosity, which 
is found to be much larger than that of the reference system, in good agreement with the Enskog results but not 
EMT. The apparent failure of EMT for the bulk viscosity may signal that the theory needs significant corrections for 
"soft" colloidal particles, or perhaps even for standard hard cores. The comparison of the MD values and theoretical 
predictions for the thermal conductivity yields the more interesting results. For heavy colloidal particles both the 
EMT and Enskog predictions appear to be in reasonable agreement with the MD simulations. On the other hand, if 
the colloidal particles are light, rric/m — 1, particularly if the size ratio is also large, the thermal conductivity is found 
to be significantly bigger than the EMT prediction. In fact, for fJc/c — 15 it appears that by adding to the solvent 
(at constant pressure) thermally insulating but rather small and light colloidal particles, the thermal conductivity is 
significantly (~ 50%) enhanced over that of the reference, pure solvent! The Enskog theory predictions also appear 
to roughly reproduce this trend. 

In conclusion, we performed MD calculations of the viscosities and thermal conductivity of a model colloidal 
suspension with colloidal particles approximately only one order of magnitude larger than the solvent molecules 
( "nanocoUoidal" suspension). We compare the calculated values with those corresponding to the pure host liquid 
("solvent") at the same pressure and temperature, and test their ratio against two conceptually distinct theories: 
effective medium theory (EMT) and Enskog transport theory. The results suggest that, quite remarkably, the standard 
EMT remains well applicable for predicting both the shear viscosity and thermal conductivity of such suspensions 
when the colloidal particles have a "typical" mass, i.e. Wc/m ~ (o'c/cr)'^. Somewhat puzzling, the available EMT 
result for the suspension bulk viscosity (essentially identical here with the reference liquid value) fails to reproduce 
the MD calculations; this may indicate that some revised theory is necessary. Estimates of the transport coefficients 
based on the Enskog transport theory are less conclusive, but appear to suggest that when applied to systems as 
the ones studied here the theory is rather inaccurate for the shear viscosity, although it may remain satisfactory for 
the bulk viscosity and the thermal conductivity. For extremely light colloidal particles, i.e. rric/m ^ 1, we find a 
significant thermal conductivity enhancement over the EMT predictions, which is roughly reproduced by the Enskog 
mixture theory. This observed mass dependence appears to distinguish in principle, from a heat conduction point of 
view, liquid suspensions from their solid phase counterparts to which EMT is most often applied. The effect may be 
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perhaps attributed to the solvent "stirring" action of the Brownian colloidal particles' motion, posited for example in 
[32j . However interesting, this behavior does not appear to be of much consequence for common, realistic suspensions, 
including nanocolloidal ones or so-called nanofluids, since for such systems the colloidal particles are dense objects 
with a mass much larger than that of the solvent molecules, and are therefore "typical" in the sense noted above. 
On the other hand the effect may be relevant for "bubbly" liquids, where the density of the host liquid is much 
larger than that of the suspended "particles" (gas bubbles). The present predictions would pertain in principle to 
nanobubbles-rich fluids, although we are not aware of any available experimental results that can either confirm or 
refute their applicability. 

I thank one of the referees for pointing out the possible relevance of the simulations with light colloidal particles to 
"bubbly" liquids. This work was performed under the auspices of the U. S. Department of Energy by University of 
California Lawrence Livermore National Laboratory under Contract No. W-7405-Eng-48. 
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TABLE I: Colloidal systems studied. 





mc/m 




4>c 


N (f) 


10 


1 


25 


0.0751 


127675 0.384 


10 


1000 


25 


0.0751 


127675 0.384 


15 


1 


20 


0.0824 


312638 0.382 


15 


1000 


20 


0.0824 


312638 0.382 


TABLE II: Comparison of the shear and bulk viscosities, and thermal conductivity (top to bottom) of the colloidal suspensions 


from MD simulations, effective medium theory (EMT) and Enskog theory for mixtures. The two EMT numbers correspond to 


volume fractions 4>c and (f 


i/c (see text). 








ac/cr 


mc/m 


MD 


EMT 


Enskog 


10 


1 


1.12 


1.19-1.25 


1.34 


10 


1000 


1.16 


1.19-1.25 


1.59 


15 


1 


1.15 


1.21-1.25 


1.54 


15 


1000 


1.21 


1.21-1.25 


1.86 


10 


1 


1.78 


1 


1.44 


10 


1000 


1.95 


1 


1.83 


15 


1 


1.86 


1 


1.71 


15 


1000 


2.20 


1 


2.22 


10 


1 


0.95 


0.85-0.89 


1.26 


10 


1000 


0.85 


0.85-0.89 


0.91 


15 


1 


1.52 


0.85-0.88 


1.40 


15 


1000 


0.80 


0.85-0.88 


0.91 
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FIG. 1: Shear viscosity autocorrelation function (normalized by its f = value) for the reference system (dashed line), and 
mixtures with cTc/c = 15; nic/m = 1000 (solid line) and mc/m = 1 (dotted line). 




FIG. 2: Autocorrelation function (normalized by its t = value) for the Lqq coefficient of the reference system (dashed line), 
and mixtures with cTc/f = 15; rric/m = 1000 (solid line) and rric/m = 1 (dotted line). 




FIG. 3: Time dependent shear viscosity (normalized by the reference system value) for mixtures with cTc/c = 15, rric/m = 1000 
(sohd line), ac/a = 10, mc/m = 1000 (dashed line), and ctc/o" = 15, mc/m = 1 (dotted line). 
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FIG. 4: Time dependent bulk viscosity (normalized by the reference system value) for mixtures with adfy = 15, malm = 1000 
(solid line), ado = 10, mc/m = 1000 (dashed line), and ac/cr = 15, mc/m = 1 (dotted line). 
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FIG. 5: Time dependent thermal conductivity (normalized by the reference system value) for mixtures with (Tc/f = 15, 
mc/m = 1000 (solid line), and ctc/o" = 10, mc/m = 1000 (dashed line). 




FIG. 6: Time dependent thermal conductivity (normalized by the reference fluid value) for the mixture with Uc/cr = 15 and 
nic/m = 1. 



